home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / jn.c < prev    next >
C/C++ Source or Header  |  1988-07-11  |  2KB  |  124 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.  The Berkeley software License Agreement
  4.  * specifies the terms and conditions for redistribution.
  5.  */
  6.  
  7. #ifndef lint
  8. static char sccsid[] = "@(#)jn.c    5.2 (Berkeley) 4/29/88";
  9. #endif /* not lint */
  10.  
  11. /*
  12.     floating point Bessel's function of
  13.     the first and second kinds and of
  14.     integer order.
  15.  
  16.     int n;
  17.     double x;
  18.     jn(n,x);
  19.  
  20.     returns the value of Jn(x) for all
  21.     integer values of n and all real values
  22.     of x.
  23.  
  24.     There are no error returns.
  25.     Calls j0, j1.
  26.  
  27.     For n=0, j0(x) is called,
  28.     for n=1, j1(x) is called,
  29.     for n<x, forward recursion us used starting
  30.     from values of j0(x) and j1(x).
  31.     for n>x, a continued fraction approximation to
  32.     j(n,x)/j(n-1,x) is evaluated and then backward
  33.     recursion is used starting from a supposed value
  34.     for j(n,x). The resulting value of j(0,x) is
  35.     compared with the actual value to correct the
  36.     supposed value of j(n,x).
  37.  
  38.     yn(n,x) is similar in all respects, except
  39.     that forward recursion is used for all
  40.     values of n>1.
  41. */
  42.  
  43. #include <math.h>
  44. #if defined(vax)||defined(tahoe)
  45. #include <errno.h>
  46. #else    /* defined(vax)||defined(tahoe) */
  47. static double zero = 0.e0;
  48. #endif    /* defined(vax)||defined(tahoe) */
  49.  
  50. double
  51. jn(n,x) int n; double x;{
  52.     int i;
  53.     double a, b, temp;
  54.     double xsq, t;
  55.     double j0(), j1();
  56.  
  57.     if(n<0){
  58.         n = -n;
  59.         x = -x;
  60.     }
  61.     if(n==0) return(j0(x));
  62.     if(n==1) return(j1(x));
  63.     if(x == 0.) return(0.);
  64.     if(n>x) goto recurs;
  65.  
  66.     a = j0(x);
  67.     b = j1(x);
  68.     for(i=1;i<n;i++){
  69.         temp = b;
  70.         b = (2.*i/x)*b - a;
  71.         a = temp;
  72.     }
  73.     return(b);
  74.  
  75. recurs:
  76.     xsq = x*x;
  77.     for(t=0,i=n+16;i>n;i--){
  78.         t = xsq/(2.*i - t);
  79.     }
  80.     t = x/(2.*n-t);
  81.  
  82.     a = t;
  83.     b = 1;
  84.     for(i=n-1;i>0;i--){
  85.         temp = b;
  86.         b = (2.*i/x)*b - a;
  87.         a = temp;
  88.     }
  89.     return(t*j0(x)/b);
  90. }
  91.  
  92. double
  93. yn(n,x) int n; double x;{
  94.     int i;
  95.     int sign;
  96.     double a, b, temp;
  97.     double y0(), y1();
  98.  
  99.     if (x <= 0) {
  100. #if defined(vax)||defined(tahoe)
  101.         extern double infnan();
  102.         return(infnan(EDOM));    /* NaN */
  103. #else    /* defined(vax)||defined(tahoe) */
  104.         return(zero/zero);    /* IEEE machines: invalid operation */
  105. #endif    /* defined(vax)||defined(tahoe) */
  106.     }
  107.     sign = 1;
  108.     if(n<0){
  109.         n = -n;
  110.         if(n%2 == 1) sign = -1;
  111.     }
  112.     if(n==0) return(y0(x));
  113.     if(n==1) return(sign*y1(x));
  114.  
  115.     a = y0(x);
  116.     b = y1(x);
  117.     for(i=1;i<n;i++){
  118.         temp = b;
  119.         b = (2.*i/x)*b - a;
  120.         a = temp;
  121.     }
  122.     return(sign*b);
  123. }
  124.